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Abstract 

A  three-dimensional  computational  fluid  dynamics  model  of  a  PEM  fuel  cell  with  serpentine  flow  field  channels  is  presented  in  this  paper. 
This  comprehensive  model  accounts  for  the  major  transport  phenomena  in  a  PEM  fuel  cell:  convective  and  diffusive  heat  and  mass  transfer, 
electrode  kinetics,  and  potential  fields.  A  unique  feature  of  the  model  is  the  implementation  of  a  voltage-to-current  (VTC)  algorithm  that 
solves  for  the  potential  fields  and  allows  for  the  computation  of  the  local  activation  overpotential.  The  coupling  of  the  local  activation 
overpotential  distribution  and  reactant  concentration  makes  it  possible  to  predict  the  local  current  density  distribution  more  accurately. 
The  simulation  results  reveal  current  distribution  patterns  that  are  significantly  different  from  those  obtained  in  studies  assuming  constant 
surface  overpotential.  Whereas  the  predicted  distributions  at  high  load  show  current  density  maxima  under  the  gas  channel  area,  low  load 
simulations  exhibit  local  current  maxima  under  the  collector  plate  land  areas. 
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1.  Introduction 

Because  of  the  highly  reactive  environment  of  a  fuel  cell  it 
is  not  possible  to  perform  detailed  in  situ  measurements  dur¬ 
ing  operation.  Yet  such  detailed  information  is  required  in 
order  to  improve  understanding  of  water  and  species  trans¬ 
port  and  of  thermal  effects.  Modeling  and  simulation  tools 
that  can  provide  such  information  and  predict  the  effect  of 
various  operating  and  material  parameters  would  contribute 
to  significantly  shorter  design  and  optimization  cycles.  The 
development  of  comprehensive  and  reliable  computational 
models  is  a  challenging  task  because  the  processes  involve 
multi-component,  multi-phase,  and  multi-dimensional  flow, 
heat  and  mass  transfer  with  electro-chemical  reactions;  all 
occurring  in  irregular  geometries  including  porous  media. 
Extensive  research  efforts  have  been  devoted  to  develop 
realistic  simulation  models  in  the  past  decade.  The  first 
two  models  were  published  in  the  early  1990s  by  Springer 
et  al.  [1]  and  Bernardi  and  Verbrugge  [2].  These  mod¬ 
els  are  one-dimensional  and  only  account  for  diffusive 
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mass  transport  and  electrochemical  kinetics.  From  1993 
to  1998,  several  two-dimensional  models  were  published 
by  Nguyen  and  White  [3],  Fuller  and  Newman  [4],  and 
Gurau  et  al.  [5].  Most  of  these  assume  some  concentration 
profile  of  reactant  species  along  the  channel  except  for 
Gurau’ s  model,  which  accounts  directly  for  convective  mass 
transport. 

The  application  of  CFD  methodology  to  fuel  cell 
modeling  has  spurred  significant  progress  including 
three-dimensional  capabilities  [6]  and  two-phase  trans¬ 
port  [7-9].  In  order  to  simplify  the  numerical  proce¬ 
dure  and  make  the  problems  computationally  tractable, 
multi-dimensional  CFD  models  have  relied  on  the  lineariza¬ 
tion  of  the  Butler- Volmer  equation  by  assuming  a  constant 
surface  overpotential.  In  addition,  most  fuel  cell  models 
presented  to  date  use  straight  gas  flow  channel  configura¬ 
tions,  whereas  serpentine  flow  channels  are  commonly  used 
in  many  fuel  cell  flow  plate  design. 

This  paper  presents  a  three-dimensional  model  with  ser¬ 
pentine  gas  flow  channels.  The  model  accounts  for  detailed 
species  mass  transport,  heat  transfer,  potential  losses  in  the 
gas  diffusion  layers  (GDF)  and  membrane,  and  electro¬ 
chemical  kinetics.  In  addition  to  the  complex  geometry,  the 
model  features  an  algorithm  that  allows  for  a  more  realistic 
representation  of  the  local  activation  overpotentials  which 
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Nomenclature 

Amea 

area  of  the  ME  A  (cm2) 

ACh 

cross  sectional  area  of  flow  channel  (cm2) 

a 

active  surface  per  unit  volume  (cm2  cm-3) 

Ch2 

local  hydrogen  concentration  (molm-3) 

/^ref 

ch2 

reference  hydrogen  concentration 
(mol  m3) 

Co2 

local  oxygen  concentration  (molm-3) 

y^ref 

co2 

reference  oxygen  concentration 
(mol  m3) 

cP 

specific  heat  capacity  (Jkg-1  K-1) 

D 

diffusion  coefficients  (m2  s-1) 

F^cell 

cell  operating  potential  (V) 

Erev 

reversible  cell  potential  (V) 

F 

Faraday’s  constant  96487  (Cmol-1) 

h 

total  specific  enthalpy  (Jkg-1) 

I 

cell  operating  current  density  (A  cm-2) 

la. 

anode  volumetric  current  density  (A  cm-3) 

k 

cathode  volumetric  current  density  (A  cm-3) 

.•ref 

lo,a 

anode  reference  exchange  current 
density  (A cm-3) 

.•ref 

Lo,c 

cathode  reference  exchange  current 
density  (Am-3) 

kp 

hydraulic  permeability  (m2) 

L 

channel  length  (m) 

Mo2 

molecular  weight  of  oxygen  (kg  mol-1) 

Mh2 

molecular  weight  of  hydrogen  (kg  mol-1) 

Mu20 

molecular  weight  of  water  (kg  mol-1) 

Vv 

net  water  flux  across  the  membrane 
(kg  m-2  s-1) 

nQ 

number  of  electrons  transfered 

P 

pressure  (Pa) 

<7 

heat  generation  (Wm-3) 

R 

universal  gas  constant  (8.314Jmol-1  K-1) 

s 

specific  entropy  (Jmol-1  K-1) 

T 

temperature  (K) 

u 

velocity  vector  (ms-1) 

Xi 

molar  fraction 

Yi 

mass  fraction 

a 

charge  transfer  coefficient 

£ 

porosity 

y 

concentration  parameter 

0 

overpotential  (V) 

A-eff 

effective  electrode  thermal 
conductivity  (Wm-1  K-1) 

Agr 

thermal  conductivity  of 
graphite  (Wm-1  K-1) 

Ae 

electrode  electronic  conductivity  (Sm-1) 

A-m 

membrane  protonic  conductivity  (Sm-1) 

Ji 

viscosity  (kgm-1  s-1) 

P 

density  (kgm-3) 

$ 

stoichiometric  flow  ratio 

leads  to  improved  prediction  of  the  local  current  density 
distribution. 


2.  Model  development 

2.7.  Model  assumptions 

In  common  with  many  detailed  three-dimensional  models, 
steady  state  operation  under  fully  humidified  conditions  is 
assumed.  While  transients  are  encountered  in  practice,  this  is 
beyond  the  scope  of  this  work.  The  model  is  also  restricted  to 
single-phase  water  transport  in  the  gas  diffusion  electrodes 
and  gas  flow  channels  and  assumes  operation  under  ideal 
heat  and  water  management  ensuring  the  membrane  remains 
fully  humidified.  Issues  related  to  two-phase  transport  and 
water  condensation  are  discussed  in  Ref.  [9]. 

Both  humidified  air  and  hydrogen  behave  as  ideal  gases 
and  since  the  characteristic  Reynolds  number  in  the  gas 
channels  are  low  the  flows  there  are  assumed  laminar. 
Other  assumptions  used  in  the  model  are  as  follows: 

•  Anode  and  cathode  gases  are  not  permitted  to  crossover, 
i.e.  the  membrane  is  impermeable. 

•  Ohmic  heating  is  neglected  as  heat  generation  is  assumed 
to  be  predominantly  associated  with  the  cathodic  electro¬ 
chemical  reaction  [10]. 

•  Water  is  assumed  to  exist  in  the  vapor  phase,  only  consis¬ 
tent  with  the  assumption  of  single-phase  water  transport. 

•  The  potential  drop  in  the  bipolar  plate  is  negligible  since 
graphite  is  a  good  conductor. 

2.2.  Computational  domain 

A  computational  model  of  an  entire  cell  would  require 
very  large  computing  resources  and  excessively  long  simula¬ 
tion  times.  The  computational  domain  in  this  study  is  there¬ 
fore  limited  to  one  turn  of  a  serpentine  channel  as  shown 
Fig.  1. 

2.3.  Governing  equations 

The  governing  Eqs.  corresponding  to  the  various  regions 
of  the  fuel  cell  are  given  below. 

The  flow  in  the  reactant  gas  channels  is  governed  by  the 
Navier-Stokes  equation. 

Continuity  equation: 

V  •  (pu)  =  0  (1) 

Momentum  equation: 

V  •  (pu  (8)  u  —  /xVu)  =  —  V(p  +  ^/xV  •  u)  +  V  •  [/x(Vu)T] 

(2) 

Energy  equation: 

V  •  ( puh )  -  V  •  (kWT)  =  0  (3) 
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Computational  Domain 


Cathode  GDL 

jTz 


Cathode  Cl 


Fig.  1.  Geometry  of  serpentine  flow  channel  fuel  cell  and  boundaries  of 
the  computational  domain.  (GDL  is  gas  diffusion  electrode;  CL  is  the 
catalyst  layer). 


The  balance  between  the  electro-osmotic  drag  of  water 
(from  anode  to  cathode)  and  back  diffusion  (from  cathode 
to  anode)  yields  the  net  water  flux  through  the  membrane 

N w  —  n(\M\{2Q  —  V  •  (pDwVTw)  (10) 

where  n&  is  the  electro-osmotic  drag  coefficient  and  Dw  is 
the  water  diffusion  coefficient  in  the  membrane. 

The  bipolar  plates  conduct  heat  and  current.  Electron  con¬ 
duction  is  assumed  to  be  very  fast  in  graphite  and  is  not 
modeled  here.  Heat  conduction  in  the  bipolar  plates  is  gov¬ 
erned  by 

V  •  (kgrV7)  =  0  (11) 

The  potential  distribution  in  the  gas  diffusion  layers  is  gov¬ 
erned  by 

V  •  (AeVV)  =  Se  (12) 


Mass  transport: 

V  •  (puYjf)  -  V  •  ( pDmWYt )  =  V  •  C oDmVYj )  (4) 

where  the  subscript  i  denotes  oxygen  at  the  cathode  side 
and  hydrogen  at  the  anode  side,  and  j  is  water  vapor  in 
both  cases.  The  diffusion  coefficients  D(a)  and  are 
Stefan-Maxwell  ternary  diffusion  coefficients  and  are  ob¬ 
tained  following  Taylor  and  Krishna  [11].  The  gas  mixture  is 
assumed  well  mixed  at  the  molecular  level,  with  all  compo¬ 
nents  sharing  the  same  velocity,  pressure,  and  temperature 
fields. 

Transport  in  the  gas  diffusion  layer  is  modeled  as  trans¬ 
port  in  a  porous  media.  The  continuity  equation  in  the  gas 
diffusion  layers  becomes: 

V  •  (pe u)  =  0  (5) 

The  momentum  equation  reduces  to  Darcy’s  law 

u  =  —  —  V  p  (6) 

and  the  mass  transport  equation  in  the  GDL  becomes 

V  •  (peuYi)  -  V  •  0 oeDuVYi )  =  V  •  ( peD^Yj )  (7) 

In  order  to  account  for  geometric  constraints  of  the  porous 
media,  the  diffusivities  are  corrected  using  the  Bruggemann 
correction  formula 

of  =  Dij  xsL5  (8) 

The  exponent  of  1.5  was  determined  empirically  [12].  It 
should  be  noted  that  the  Bruggemann  correction  and  the  al¬ 
ternative  correction  based  on  tortuosity  are  essentially  equiv¬ 
alent  for  low  tortuosities  and  porosities  in  the  range  0.4-0. 5. 
Heat  transfer  in  the  in  the  gas  diffusion  layers  is  governed 

by, 

V  •  (peuh)  -  V  •  (fceff£V7)  =  e  Sq  (9) 

where  Sq  is  the  source  term  for  heat  exchange  to  and  from 
the  solid  matrix  of  the  porous  media. 


Potential  loss  in  the  membrane  is  due  to  resistance  to 
proton  transport  across  the  membrane,  and  is  governed  by  a 
similar  equation, 

v  .  (AmVV)  =  (13) 

2.4.  Electrochemical  reaction  kinetics 


The  consumption  of  reactant  species  and  the  production 
of  water  and  heat  are  modeled  as  sink  and  source  terms  in 
the  catalyst  layers.  The  mass  consumption  rate  of  oxygen 
per  unit  volume  is  given  by: 


(14) 


The  production  of  water  is  modeled  as  a  source  term  based 
on  the  local  current  density, 


(15) 


At  the  anode  catalyst  layer,  hydrogen  is  consumed  to  pro¬ 
duce  electrons  and  protons.  The  consumption  of  hydrogen 
is  given  by 


Svh  =  - 


Mh2  . 

- -i- 

2  F 


(16) 


In  this  model,  heat  generation  is  assumed  to  be  predomi¬ 
nantly  due  to  the  electrochemical  reactions,  and  ohmic  heat¬ 
ing  is  not  currently  accounted  for.  Furthermore,  Lampinen 
and  Fomino  [10]  showed  that  heat  generation  from  the  an¬ 
ode  reaction  is  negligible  compared  to  the  cathode  reaction, 
and  hence  only  cathodic  heat  generation  is  considered: 


T(-As) 

ntF 


+  7?  act 


ic 


(17) 


where  T  is  the  local  temperature,  As  is  the  entropy  of  the 
chemical  reaction,  nQ  is  the  number  of  electrons  transferred 
per  mole  of  hydrogen,  r]act  is  the  activation  overpotential. 
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The  local  current  density  distribution  in  the  catalyst  layers 
is  modeled  by  the  Butler- Volmer  equation. 


i 


c 


•ref  /  Cq2 

0,c  \  f 

\co2 


ro2 


exp 


—  exp 


•  _  -ref  [  ^H2 
a  o,a  l  ^ref 

\CH2 


Kh2 


exp 


—  exp 


(18) 


(19) 


where  iT*fc  and  are  the  cathode  reference  exchange  cur¬ 
rent  density;  and  ac  are  the  cathode  anodic  and  cathodic 
charge  transfer  coefficients  respectively;  and  y  is  an  empiri¬ 
cal  concentration  parameter.  The  values  of  the  electrochem¬ 
ical  transport  parameters  are  taken  from  [13]  and  are  listed 
in  Table  1.  The  variation  of  the  exchange  current  density 
with  temperature  was  computed  using  the  empirical  relation 
given  by  Parthasarathy  et  al.  [14]. 

Given  an  input  cell  potential  Ece n,  the  local  cathode  acti¬ 
vation  overpotential  is  calculated  from 


7?act,c  —  ^rev  ^cell  ^ohmic,gdl  protonic,  mem  ^?act,a 


(20) 


where  Erev  is  obtained  from  the  Nernst  equation 


EKY  =  E0-^ln(^^\  (21) 

2  F  W AJ 

The  ohmic  and  protonic  overpotentials  in  Eq.  (20)  are 
calculated  from  the  potential  Eqs.  (12)  and  (13)  using 

>7  ohmic  =  VCL  Vref  (22) 

Volume  —  Ecl,c  —  P(2L,a  (23) 

where  Vcl  is  the  potential  at  the  catalyst  layer  and  Vref  is 
the  reference  potential  taken  at  the  interface  between  the  gas 
diffusion  layer  and  the  bipolar  plate. 

The  anode  activation  overpotential  is  small  and  assumed 
to  be  uniform.  It  is  calculated  from  a  conditional  loop.  An 
initial  guess  of  the  anode  activation  overpotential  is  obtained 
from  the  average  local  current  density  at  the  anode  side,  /a 
using  the  Butler- Volmer  equation  (Eq.  (19)).  An  updated 
value  of  the  anode  activation  is  then  computed  based  on  the 
charge  conservation  condition 

Nc  Na 

^total  =  TjicJ  X  Vj  =  X  Vj  (24) 

7=1  7=1 

where  Nc  and  Va  are  the  total  number  of  computational  cells 
in  the  cathode  and  anode  catalyst  layers  respectively;  Vj  is 
the  computational  cell  volume. 

It  should  be  noted  that  in  many  previous  studies  the  cur¬ 
rent  is  prescribed  and  a  constant  surface  overpotential  is 
assumed  in  order  to  linearize  the  Butler- Volmer  equation. 


Table  1 


Model  parameters 


Parameter 

Symbol 

Value 

Channel  length  (mm) 

L 

50 

Channel  width  (mm) 

w 

1 

Channel  height  (mm) 

h 

1 

Land  area  width  (mm) 

W\ 

1 

Gas  diffusion  layer  thickness  (pm) 

Cell 

260 

Wet  membrane  thickness  (Nation  117)  (pm) 

t  mem 

230 

Catalyst  layer  thickness  (pm) 

Cl 

28 

Membrane  porosity 

£ 

0.28 

Electrode  porosity 

£ 

0.4 

Membrane  ionic  conductivity  (humidified  Nation®  117)  (Sm-1) 

Cm 

17 

Electrode  electronic  conductivity  (estimate)  (Sm-1) 

Ce 

570 

Cathodic  charge  transfer  coefficient  for  cathode  reaction 

01 c,c 

2 

Anodic  charge  transfer  coefficient  for  cathode  reaction 

0i  a,c 

2 

Cathodic  charge  transfer  coefficient  for  anode  reaction 

OL  c,a 

1 

Anodic  charge  transfer  coefficient  for  anode  reaction 

aa,  a 

1 

Hydrogen  reference  concentration  (mol cm-3) 

y^ref 

H2 

5.64  E 

-5 

Oxygen  reference  concentration  (mol  cm-3) 

r'xf 

c02 

3.39  E 

-6 

Oxygen  concentration  parameter 

Ko2 

0.5 

Hydrogen  concentration  parameter 

Kh2 

0.25 

Cathode  reference  exchange  current  density  (A cm-3) 

7-ref 

lo,c 

1.0  E- 

5 

Anode  reference  exchange  current  density  (A cm-3) 

?-ref 

G,a 

1.4  E5 

Membrane  thermal  conductivity  (W  m- 1  K- 1 ) 

k] m 

0.455 

Electrode  thermal  conductivity  (Ballard  AvCarb®-P150)  (Wm-1  K-1) 

Ceff 

1.3 

Membrane  hydraulic  permeability  (m2) 

kp 

1.8  E- 

18 

Electrode  hydraulic  permeability  (m2) 

kp 

4.73  E- 

-19 
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Table  2 


Base  case  operating  conditions 


Parameter 

Symbol 

Value 

Anode  pressure 

Pa 

3  atm 

Cathode  pressure 

Pc 

3  atm 

Cell  temperature 

T 

353  K 

Relative  humidity  of  inlet  fuel 

(Pa 

100% 

Relative  humidity  of  inlet  air 

(Pc 

100% 

Air  stoichiometric  ratio 

$ c 

2 

Fuel  stoichiometric  ratio 

2 

In  this  work,  the  local  activation  overpotential  distribution 
is  calculated  by  using  the  iterative  Voltage-to-Current  algo¬ 
rithm  described  in  Section  3. 


2.5.  Boundary  conditions 


The  inlet  velocities  of  air  and  fuel  are  calculated  based 
on  the  cell  current  and  stoichiometric  flow  ratio  as  follow, 


U inlet 


-A  ME  A  RT 

neF  ACh  xiPi 


(25) 


where  nQ  equals  2  for  the  anode  flow  side  and  4  for  the 
cathode  side. 

Pressure  boundary  conditions  are  prescribed  at  the  outlets 
to  simulate  the  operating  pressure  of  the  fuel  cell.  Symme¬ 
try  boundary  conditions  are  applied  in  the  x-y  plane  of  the 
computational  domain.  Zero  heat  flux  is  applied  at  the  x-z 
plane  of  the  conducting  boundary  surfaces. 

A  combination  of  Dirichlet  and  Neumann  boundary  con¬ 
ditions  are  used  to  solve  the  electronic  and  protonic  potential 


equations.  A  reference  potential  is  applied  at  the  interface 
between  the  GDLs  and  the  bipolar  plate  (land  area). 

The  geometric  and  physico-chemical  parameters  are  listed 
in  Table  1,  and  the  base  case  operating  conditions  in  Table  2. 


3.  Solution  algorithm 

The  governing  equations  and  appropriate  boundary  con¬ 
ditions  were  implemented  and  solved  using  a  3  D  commer¬ 
cial  Computational  Fluid  Dynamics  package  (CFX-4.3)  that 
employs  a  finite  volume  formulation.  The  discretized  sets 
of  algebraic  equations  are  solved  using  the  full  field  Stone’s 
method  for  the  momentum  and  scalar  transport  equations, 
and  an  algebraic  multi-grid  method  for  the  pressure  and  en¬ 
ergy  equations.  The  properties  are  updated  after  each  global 
iterative  loop  based  on  the  new  local  gas  composition  and 
temperature  using  the  relations  given  in  Ref.  [6]. 

The  computational  mesh  consisted  of  a  body-fitted  grid 
with  351000  computational  cells.  The  Fortran  user  subrou¬ 
tines  developed  previously  [6]  to  solve  the  electrochemical 
component  of  the  model  were  adapted  and  a  new  VTC  al¬ 
gorithm  was  implemented. 

A  unique  feature  of  the  iterative  VTC  algorithm  devel¬ 
oped  in  this  work  is  its  capability  for  accurate  calculation  of 
the  local  activation  overpotentials,  which  in  turn  results  in 
improved  prediction  of  the  current  density  distribution.  This 
algorithm  is  capable  of  predicting  the  cell  current  density 
based  on  a  target  cell  operating  voltage.  The  flowchart  for 
the  algorithm  is  shown  in  Fig.  2. 

The  improved  predictive  capabilities  of  the  iterative  VTC 
algorithm  come  with  the  penalty  of  higher  computational 


Fig.  2.  Voltage-to-current  (VTC)  algorithm. 
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Fig.  3.  Comparison  of  model  and  experimental  polarization  curves. 


costs.  Depending  on  the  prescribed  operating  voltage,  sim¬ 
ulations  with  the  VTC  algorithm  require  about  6000-8000 
iterations  to  achieve  convergence;  the  higher  the  load  the 
slower  the  convergence. 


ETA 


© 

o 


-0.03  -0.02  -0.01 

Cell  Length  (m) 


0.3724 

0.3659 

0.3593 

0.3528 

0.3463 

0.3398 

0.3332 

0.3267 

0.3202 

0.3137 


4.  Results  and  discussions 

The  simulation  results  for  base  case  operating  conditions 
were  verified  against  measurements  of  Ref.  [15]  obtained  for 
the  same  conditions  (Table  2).  The  computed  polarization 
curve  shown  in  Fig.  3  is  in  good  agreement  with  the  exper¬ 
imental  polarization  curve  in  the  low  load  region.  However, 
the  model  cell  current  densities  in  the  mass  transport  limited 
region  (>1.25  A  cm-2)  are  higher  than  the  experimental  val¬ 
ues.  This  discrepancy  is  a  common  feature  of  single-phase 
models  where  the  effect  of  reduced  oxygen  transport  due 
to  water  flooding  at  the  cathode  at  high  current  density  [9] 


Fig.  5.  Activation  overpotential  (in  volts)  distribution  in  the  first  catalyst 
layer  for  two  loading  conditions:  (a)  0.3  Acm~2  and  (b)  1.2  A  cm-2. 


cannot  be  accounted  for.  In  addition  to  this  flooding  effect, 
anode  drying  can  also  be  a  contributing  factor  to  the  reduced 
performance  at  high  current  density. 

The  detailed  distribution  of  oxygen  mass  fractions  for  two 
different  loading  conditions  is  shown  in  Fig.  4.  Oxygen  con¬ 
centration  decreases  from  the  inflow  channel  to  the  outflow 
channel.  In  the  GDL,  oxygen  concentration  under  the  land 
area  is  smaller  than  that  under  the  channel  area.  This  effect 
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Fig.  4.  Oxygen  molar  fraction  distribution  in  the  cathode  side  for  two  loading  conditions:  (a)  0.3  A  cm  2  and  (b)  1.2  A  cm 
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is  more  pronounced  for  high  load  conditions,  where  oxygen 
under  the  land  area  is  almost  completely  depleted. 

The  variation  of  the  cathode  activation  overpotentials  is 
shown  in  Fig.  5.  For  both  loading  conditions,  the  distribution 
patterns  of  activation  overpotentials  are  similar,  with  higher 
values  under  the  land  area.  This  uneven  distribution  is  due 
to  lower  ohmic  losses  under  the  land  area  than  under  the 
channel  area.  Since  the  activation  overpotential  has  an  ex¬ 
ponential  effect  on  the  magnitude  of  local  current,  the  con¬ 
ditions  for  generating  current  under  the  land  area  are  more 
favorable  in  the  absence  of  oxygen  transport  limitations.  The 
enhanced  capability  for  producing  current  under  the  land 
areas  when  oxygen  transport  is  not  limited  is  illustrated  in 
Fig.  6  which  shows  distribution  patterns  that  are  quite  dif¬ 
ferent  from  those  obtained  in  previous  studies  in  which  a 
constant  overpotential  is  assumed.  For  the  low  load  case  il¬ 
lustrated,  oxygen  diffusion  is  not  limiting,  and  the  current 
density  is  higher  under  the  land  area.  However,  as  loading 
increases,  oxygen  transport  limitations  under  the  land  area 
become  significant,  resulting  in  the  shift  of  higher  current 
density  towards  the  channel  area  where  oxygen  concentra¬ 
tion  is  higher. 

Ohmic  overpotential  is  the  loss  associated  with  resistance 
to  electron  transport  in  the  gas  diffusion  layers.  For  a  given 
load,  the  magnitude  of  this  overpotential  is  dependent  on 
the  path  of  the  electrons.  The  longer  the  path,  the  larger  the 
potential  drop.  Under  the  channel  area  the  catalyst  layer  is 
farthest  from  the  collector  plate  and  locally  ‘edge’  collection 
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Fig.  6.  Relative  current  density  distribution  in  the  first  catalyst  layer 
for  two  loading  conditions:  (a)  0.3  A  cm-2  and  (b)  1.2  A  cm-2  (current 
densities  normalized  by  100  A  cm-3). 


takes  place,  i.e.  the  current  flows  along  a  longer  path  along 
the  width  (z-direction)  rather  than  across  the  GDL.  Ohmic 
losses  are  thus  higher  than  under  the  gas  diffusion  layer. 
Fig.  7  shows  ohmic  overpotential  distribution  patterns  that 
are  similar  for  both  loading  conditions.  However,  the  mag¬ 
nitude  of  the  potential  loss  increases  with  cell  loading. 

As  indicated  earlier,  heat  generation  is  only  considered 
from  the  electrochemical  reaction  on  the  cathode  side  since 
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Fig.  7.  Ohmic  overpotential  in  the  cathode  gas  diffusion  layer  and  channel  for  two  loading  conditions:  (a)  0.3  A  cm  2  (b)  1.2  A  cm 
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the  activation  is  negligible  in  the  anode  side.  The  temperature 
distribution  inside  a  fuel  cell  is  highly  dependent  on  the 
loading  conditions  and  associated  electrochemical  activity. 
Temperature  distributions  for  the  intermediate  and  high  load 
conditions  are  shown  in  Fig.  8.  In  both  the  cases,  the  highest 
temperatures  are  located  at  the  cathode  catalyst  layer.  Since 
the  membrane  conductivity  is  quite  low  and  there  is  no  heat 
generation  on  the  anode  side,  the  anode  side  temperature  is 
quite  uniform  and  equal  to  the  anode  gas  stream  temperature. 
On  the  other  hand,  the  cathode  gas  temperature  close  to  the 
channel  gas  diffusion  layer  interface  is  from  1  to  7  °C  higher 
than  the  nominal  operating  temperature. 

5.  Conclusions 

A  three-dimensional  computational  fluid  dynamics  model 
of  a  PEM  fuel  cell  with  serpentine  flow  channels  was  devel¬ 
oped.  This  model  provides  valuable  information  about  the 
transport  phenomena  inside  the  fuel  cell  such  as  reactant  gas 
concentration  distribution,  temperature  distribution,  poten¬ 
tial  distribution  in  the  membrane  and  gas  diffusion  layers, 
activation  overpotential  distribution,  and  local  current  den¬ 
sity  distribution. 

A  unique  feature  of  this  model  is  the  implementation  of 
the  voltage-to-current  algorithm  that  allows  for  a  more  re¬ 
alistic  spatial  variation  of  the  electrochemical  kinetics.  In 
addition,  the  three-dimensional  activity  of  the  catalyst  layer 
is  also  accounted  for  in  this  model. 

The  computational  procedure  involves  the  coupling  of 
the  potential  field  with  the  reactant  species  concentration 


field.  Rather  than  assuming  a  constant  surface  overpotential 
over  the  entire  fuel  cell,  the  spatial  variation  of  the  cathodic 
surface  overpotential  is  computed  locally,  resulting  in  im¬ 
proved  prediction  of  the  local  current  density  distribution. 
The  current  density  distribution  patterns  are  found  to  vary 
with  loading  conditions.  At  low  load,  the  current  density 
is  higher  under  the  land  area.  As  the  load  is  increased,  the 
current  density  maxima  shift  towards  the  center  of  the  chan¬ 
nel.  These  local  current  density  distribution  patterns  are 
radically  different  from  those  obtained  with  models  that  do 
not  account  for  the  non-uniformity  of  surface  overpotential. 
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